function [CRL,CRR,CRB,CRT,totLL,totLR,totLB,totLT] = ...
    cresp_mc(phi,dx,dy,dc,mt,xfm,yfm,numg,alphaB,alL,alR,alT)
% Current response generator for mesh center.  The flux is an input,
% and we compute the outgoing partial currents for transmission to other
% cells and the net current for leakage.

% number of fine meshes in x and y coordinates
N = sum(xfm);
M = sum(yfm);

if nargin==9
    alL = zeros(M,numg);
    alR = alL; alT = 0*alphaB;
end

% number of coarse meshes in x and y coordinates
CN = length(xfm);
CM = length(yfm); 
LL = zeros(M,1); LR = LL; LB = zeros(N,1); LT = LB;

% coarse mesh index for fine mesh
cix=zeros(N,1);
ciy=zeros(M,1);
j = 0;
for i = 1:CN
    cix( (j+1):(j+xfm(i)) )=i;
    j = sum(xfm(1:i));
end
j = 0;
for i = 1:CM
    ciy( (j+1):(j+yfm(i)) )=i;
    j = sum(yfm(1:i));
end

CRL = zeros(M,numg); CRR = CRL;
CRB = zeros(N,numg); CRT = CRB;

for g = 1:numg
    
    % for L and R
    for j = 1:M

        % compute left response:  Jout = 0.25*phi + 0.5*D*(phi'-phi)/dx
        i   = 1;
        k   = i+(j-1)*( N );
        mC  = mt( cix(i), ciy(j) );  
        %tmp =  (4*dc(mC,g)-dx(1))/(4*dc(mC,g)+dx(i));
        %CRL( j, g) = 0.25*phi(k,g)*(1+tmp);
        %LL(j)      = LL(j) + 2*dc(mC,g)*(phi(k,g))/(4*dc(mC,g)+dx(i));
        phi_prime_old   = 2*phi(k,g)/(4*dc(mC,g)+dx(i))-8*alL(j,g)/(4*dc(mC,g)+dx(i));
        LL(j)           = LL(j) + dc(mC,g)*phi_prime_old;     
        CRL( j, g)      = 0.25*phi(k,g) + (0.5*dc(mC,g)-0.125*dx(1))*phi_prime_old;
   
%         % f-d approach        
%         mA = mt(cix(i+1),ciy(j)); mB = mt(cix(i+1),ciy(j));
%         kA = (i+1)+(j-1)*( N );   kB = (i+2)+(j-1)*( N );
%         invdx     = 1/dx(1);
%         phi_prime = invdx*(-2*phi(k,g)+3*phi(kA,g)-phi(kB,g));
%         LL(j)     = LL(j) + dc(mC,g)*phi_prime;
%         CRL(j,g)  = 0.25*phi(k,g) + (0.5*dc(mC,g)-0.125*dx(1))*phi_prime;
%         
        % compute right response
        i = N;
        k = i+(j-1)*( N );
        mC = mt( cix(i), ciy(j) );  
        tmp =  (4*dc(mC,g)-dx(i))/(4*dc(mC,g)+dx(i));
        phi_prime_old = -2*phi(k,g)/(4*dc(mC,g)+dx(i))+ ...
                        8*alR(j,g)/(4*dc(mC,g)+dx(i));
        %LR(j)      = LR(j) + 2*dc(mC,g)*(phi(k,g))/(4*dc(mC,g)+dx(i));          
        CRR( j, g) = 0.25*phi(k,g)*(1+tmp);
        LR(j)      = LR(j) - dc(mC,g)*phi_prime_old;          
        CRR( j, g) = 0.25*phi(k,g) - (0.5*dc(mC,g)-0.125*dx(1))*phi_prime_old;
    end
    
    % for B and T
    for i = 1:N

        % compute bottom response
        j = 1;
        k = i+(j-1)*( N );
        mC         = mt( cix(i), ciy(j) );  
        %if 1==0
        phi_prime_old = 2*phi(k,g)/(4*dc(mC,g)+dy(j)) - 8*alphaB(i,g)/(4*dc(mC,g)+dy(j));
%         tmp        =  (4*dc(mC,g)-dy(j))/(4*dc(mC,g)+dy(j));
%         CRB( i, g) = 0.25*phi(k,g)*(1+tmp) - 2*alphaB(j,g)*tmp;
%         LB(i)      = LB(i) + 2*dc(mC,g)*(phi(k,g)-2*alphaB(j,g))/(4*dc(mC,g)+dy(j));    
        LB(i)       = LB(i) + dc(mC,g)*phi_prime_old;  
        CRB(i,g)    = 0.25*phi(k,g) + (0.5*dc(mC,g)-0.125*dy(1))*phi_prime_old;
        %else
%         mA = mt(cix(i),ciy(j+1)); mB = mt(cix(i),ciy(j+2));
%         kA = (i)+(j-1+1)*( N );   kB = (i)+(j-1+2)*( N );
%         invdy     = 1/dy(1);
%         phi_prime = invdy*(-2*phi(k,g)+3*phi(kA,g)-phi(kB,g));
%         LB(i)     = LB(i) + dc(mC,g)*phi_prime;
%         CRB(i,g)  = 0.25*phi(k,g) + (0.5*dc(mC,g)-0.125*dy(1))*phi_prime;
        %end
        % compute top response
        j = M;
        k = i+(j-1)*( N ); 
        mC         = mt( cix(i), ciy(j) );  
        tmp        = (4*dc(mC,g)-dy(j))/(4*dc(mC,g)+dy(j));
        phi_prime_old = -2*phi(k,g)/(4*dc(mC,g)+dy(j)) + ...
                        8*alT(i,g)/(4*dc(mC,g)+dy(j));
        %CRT( i, g) = 0.25*phi(k,g)*(1+tmp);
        %LT(i)      = LT(i) + 2*dc(mC,g)*(phi(k,g))/(4*dc(mC,g)+dy(j));
        LT(i)       = LT(i) - dc(mC,g)*phi_prime_old;  
        CRT(i,g) 	= 0.25*phi(k,g) - (0.5*dc(mC,g)-0.125*dy(1))*phi_prime_old;
    end
    
end
%dxx = [0.5*dx(1) dx(1:end-1)' 0.5*dx(end)]';
%dyy = [0.5*dy(1) dy(1:end-1)' 0.5*dy(end)]';
totLL = sum( dy.*( LL ) );
totLR = sum( dy.*( LR ) );
totLB = sum( dx.*( LB ) );
totLT = sum( dx.*( LT ) );

poop=1;
end
